# Replication code for Daniel Lobo and Ryan Brutger's  manuscript (Last updated October 20, 2024)
# "Fairness According to Whom"
# Code for the appendix is in the companion R script

# Code was run on R version 4.0.1 (2020-06-06) -- "See Things Now"
# Using MAcBook Pro using macOS Catalina, with 2.9GHz Dual-Core Intel Core i5

# Race/Ethnicity Codings
# 1= White
# 2= African American
# 3 = Hispanic
# 4 = Asian
# 5 = Native American
# 6 = Pacific Islander
# 7 = Other

################# 1. Load libraries, load and clean data ##############

# First, set your working directory using the setwd() command:
setwd()

data <- get(load("Fairness_to_Whom_Replication.RData"))

#Load required R packages (download first if necessary)
library(foreign)
library(Hmisc)
library(ggplot2)
library(psych)
library(stargazer)
library(mediation)

# Create indicator for white respondents
data$white <- ifelse(data$race==1, 1, 0)
# create indicator for men
data$men <- ifelse(data$gender==1, 1, 0)

# Subset to only white and Black respondents
data <- subset(data[data$race<3, ])


## Analyze the Trade Concessions Study
## "V1" is used to indicate the data is from the Trade Concessions study

# Create data subset to those who answered the support dependent variable
dataV1 <- subset(data, is.na(data$V1_supp)==FALSE ) 
# Create data subset those who answered fairness question
dataV1_fair <- subset(dataV1[is.na(dataV1$V1_fair)==FALSE, ])

# Rescaling function for more easily interpretable results
rescale <- function(x){
  return((x-min(x,na.rm=TRUE))/(max(x-min(x,na.rm=TRUE),na.rm=TRUE)))
}

# Scale cooperative internationalism and  national attachment measures for easy interpretability
dataV1_fair$ci <- rescale(dataV1_fair$ci)
dataV1_fair$NatAtt <- rescale(dataV1_fair$NatAtt)


# Set critical value used for confidence intervals
c.value <-qnorm(.95)

# Create measure of those expressing positive fairness
dataV1_fair$dich_fair <- ifelse(dataV1_fair$V1_fair>0, 1, 0)

# Subset data by respondent's race
dataV1W_fair <- subset(dataV1_fair, dataV1_fair$race==1 ) 
dataV1B_fair <- subset(dataV1_fair, dataV1_fair$race==2 ) 

# Table 1 of manuscript (presentation of variables is re-ordered in manuscript)
mod.1 <- lm(V1_fair ~ V1_equal*white + V1_fav*white, data=dataV1_fair)
stargazer(mod.1)
summary(mod.1) 
# Outputs from summary model are reported in the text of manuscript

# Create Figure for white and Black's fairness scores by treatment
# Figure 1 of manuscript: export to PDF as 7x4

# Generate estimates for white respondents
V1_equal <- mean(dataV1W_fair$V1_fair[dataV1W_fair$V1_equal==1])
se.V1_equal <- sqrt(var(dataV1W_fair$V1_fair[dataV1W_fair$V1_equal==1])/length(dataV1W_fair$V1_fair[dataV1W_fair$V1_equal==1])) 
ci.up.V1_equal <- V1_equal + se.V1_equal * c.value   
ci.low.V1_equal <- V1_equal  - se.V1_equal * c.value 

V1_unfav <- mean(dataV1W_fair$V1_fair[dataV1W_fair$V1_unfav==1])
se.V1_unfav <- sqrt(var(dataV1W_fair$V1_fair[dataV1W_fair$V1_unfav==1])/length(dataV1W_fair$V1_fair[dataV1W_fair$V1_unfav==1])) 
ci.up.V1_unfav <- V1_unfav + se.V1_unfav * c.value   
ci.low.V1_unfav <- V1_unfav  - se.V1_unfav * c.value

V1_fav <- mean(dataV1W_fair$V1_fair[dataV1W_fair$V1_fav==1])
se.V1_fav <- sqrt(var(dataV1W_fair$V1_fair[dataV1W_fair$V1_fav==1])/length(dataV1W_fair$V1_fair[dataV1W_fair$V1_fav==1])) 
ci.up.V1_fav <- V1_fav + se.V1_fav * c.value   
ci.low.V1_fav <- V1_fav  - se.V1_fav * c.value

# Consolidate fairness estimates for plot
V1.means.Wf <- c(V1_equal, V1_unfav, V1_fav)
V1.ci.up.Wf <- c(ci.up.V1_equal, ci.up.V1_unfav, ci.up.V1_fav)
V1.ci.low.Wf <- c(ci.low.V1_equal, ci.low.V1_unfav, ci.low.V1_fav)
V1.names.Wf <- c("Equal", "Unfavorable", "Favorable")

# Generate estimates for Black respondents
V1_equal <- mean(dataV1B_fair$V1_fair[dataV1B_fair$V1_equal==1])
se.V1_equal <- sqrt(var(dataV1B_fair$V1_fair[dataV1B_fair$V1_equal==1])/length(dataV1B_fair$V1_fair[dataV1B_fair$V1_equal==1])) 
ci.up.V1_equal <- V1_equal + se.V1_equal * c.value   
ci.low.V1_equal <- V1_equal  - se.V1_equal * c.value 

V1_unfav <- mean(dataV1B_fair$V1_fair[dataV1B_fair$V1_unfav==1])
se.V1_unfav <- sqrt(var(dataV1B_fair$V1_fair[dataV1B_fair$V1_unfav==1])/length(dataV1B_fair$V1_fair[dataV1B_fair$V1_unfav==1])) 
ci.up.V1_unfav <- V1_unfav + se.V1_unfav * c.value   
ci.low.V1_unfav <- V1_unfav  - se.V1_unfav * c.value

V1_fav <- mean(dataV1B_fair$V1_fair[dataV1B_fair$V1_fav==1])
se.V1_fav <- sqrt(var(dataV1B_fair$V1_fair[dataV1B_fair$V1_fav==1])/length(dataV1B_fair$V1_fair[dataV1B_fair$V1_fav==1])) 
ci.up.V1_fav <- V1_fav + se.V1_fav * c.value   
ci.low.V1_fav <- V1_fav  - se.V1_fav * c.value

# Consolidate fairness estimates for plot
V1.means.Bf <- c(V1_equal, V1_unfav, V1_fav)
V1.ci.up.Bf <- c(ci.up.V1_equal, ci.up.V1_unfav, ci.up.V1_fav)
V1.ci.low.Bf <- c(ci.low.V1_equal, ci.low.V1_unfav, ci.low.V1_fav)
V1.names.Bf <- c("Equal", "Unfavorable", "Favorable")

# Figure 1 of manuscript: White/Black Fairness scores by Equal, Unfav, Fav treatments
# Export to PDF as 7x4
par(mar=c(3,6,1,2))
plot(1:3, V1.means.Wf , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.35, .9), xlim = c(0.75, 3.25), pch = 22, type = "p", xaxt="n", cex=1)
mtext("Average Fairness ", side=2, at=.25, cex=1.9, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(V1.means.Wf)){
  lines(c(i,i), c(V1.ci.up.Wf[i], V1.ci.low.Wf[i]))
  mtext(V1.names.Wf[i], side=1, at=i, 1.12, cex=1.9)
}
points(1.15:3.15, V1.means.Bf, pch=19, cex=1.5)
for(i in 1:length(V1.means.Bf)){
  lines(c(i+.15,i+.15), c(V1.ci.up.Bf[i], V1.ci.low.Bf[i]))
}
text(3.01, .85, labels = "Whites", cex=1.5)
text(3, .75, labels = "Blacks", cex=1.5)
points(2.7, .85, pch=22, cex=1.5)
points(2.7, .75, pch=19, cex=1.5)

# Create Figure for white and Black respondents' support scores by treatment
# Figure 2 of manuscript: export to PDF as 7x4
dataV1W_supp <- subset(dataV1, dataV1$race==1 ) 
dataV1B_supp <- subset(dataV1, dataV1$race==2 ) 

# Generate estimates for white respondents
V1_equal <- mean(dataV1W_supp$V1_supp[dataV1W_supp$V1_equal==1])
se.V1_equal <- sqrt(var(dataV1W_supp$V1_supp[dataV1W_supp$V1_equal==1])/length(dataV1W_supp$V1_supp[dataV1W_supp$V1_equal==1])) 
ci.up.V1_equal <- V1_equal + se.V1_equal * c.value   
ci.low.V1_equal <- V1_equal  - se.V1_equal * c.value 

V1_unfav <- mean(dataV1W_supp$V1_supp[dataV1W_supp$V1_unfav==1])
se.V1_unfav <- sqrt(var(dataV1W_supp$V1_supp[dataV1W_supp$V1_unfav==1])/length(dataV1W_supp$V1_supp[dataV1W_supp$V1_unfav==1])) 
ci.up.V1_unfav <- V1_unfav + se.V1_unfav * c.value   
ci.low.V1_unfav <- V1_unfav  - se.V1_unfav * c.value

V1_fav <- mean(dataV1W_supp$V1_supp[dataV1W_supp$V1_fav==1])
se.V1_fav <- sqrt(var(dataV1W_supp$V1_supp[dataV1W_supp$V1_fav==1])/length(dataV1W_supp$V1_supp[dataV1W_supp$V1_fav==1])) 
ci.up.V1_fav <- V1_fav + se.V1_fav * c.value   
ci.low.V1_fav <- V1_fav  - se.V1_fav * c.value

# Consolidate support estimates for plot
V1.means.Ws <- c(V1_equal, V1_unfav, V1_fav)
V1.ci.up.Ws <- c(ci.up.V1_equal, ci.up.V1_unfav, ci.up.V1_fav)
V1.ci.low.Ws <- c(ci.low.V1_equal, ci.low.V1_unfav, ci.low.V1_fav)
V1.names.Ws <- c("Equal", "Unfavorable", "Favorable")

# Generate estimates for Black respondents
V1_equal <- mean(dataV1B_supp$V1_supp[dataV1B_supp$V1_equal==1])
se.V1_equal <- sqrt(var(dataV1B_supp$V1_supp[dataV1B_supp$V1_equal==1])/length(dataV1B_supp$V1_supp[dataV1B_supp$V1_equal==1])) 
ci.up.V1_equal <- V1_equal + se.V1_equal * c.value   
ci.low.V1_equal <- V1_equal  - se.V1_equal * c.value 

V1_unfav <- mean(dataV1B_supp$V1_supp[dataV1B_supp$V1_unfav==1])
se.V1_unfav <- sqrt(var(dataV1B_supp$V1_supp[dataV1B_supp$V1_unfav==1])/length(dataV1B_supp$V1_supp[dataV1B_supp$V1_unfav==1])) 
ci.up.V1_unfav <- V1_unfav + se.V1_unfav * c.value   
ci.low.V1_unfav <- V1_unfav  - se.V1_unfav * c.value

V1_fav <- mean(dataV1B_supp$V1_supp[dataV1B_supp$V1_fav==1])
se.V1_fav <- sqrt(var(dataV1B_supp$V1_supp[dataV1B_supp$V1_fav==1])/length(dataV1B_supp$V1_supp[dataV1B_supp$V1_fav==1])) 
ci.up.V1_fav <- V1_fav + se.V1_fav * c.value   
ci.low.V1_fav <- V1_fav  - se.V1_fav * c.value

# Consolidate support estimates for plot
V1.means.Bs <- c(V1_equal, V1_unfav, V1_fav)
V1.ci.up.Bs <- c(ci.up.V1_equal, ci.up.V1_unfav, ci.up.V1_fav)
V1.ci.low.Bs <- c(ci.low.V1_equal, ci.low.V1_unfav, ci.low.V1_fav)
V1.names.Bs <- c("Equal", "Unfavorable", "Favorable")

# Figure 2 of manuscript: export to PDF as 7x4
par(mar=c(3,6,1,2))
plot(1:3, V1.means.Ws , main = "", 
     ylab = "", xlab = "", cex.lab=2,
     ylim = c(-.35, .9), xlim = c(0.75, 3.25), pch = 22, type = "p", xaxt="n", cex=1)
mtext("Average Support ", side=2, at=.25, cex=1.9, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(V1.means.Ws)){
  lines(c(i,i), c(V1.ci.up.Ws[i], V1.ci.low.Ws[i]))
  mtext(V1.names.Ws[i], side=1, at=i, 1.12, cex=1.9)
}
points(1.15:3.15, V1.means.Bs, pch=19, cex=1.5)
for(i in 1:length(V1.means.Bs)){
  lines(c(i+.15,i+.15), c(V1.ci.up.Bs[i], V1.ci.low.Bs[i]))
}
text(3.01, .85, labels = "Whites", cex=1.5)
text(3, .75, labels = "Blacks", cex=1.5)
points(2.7, .85, pch=22, cex=1.5)
points(2.7, .75, pch=19, cex=1.5)

# Calculate support quantities mentioned in text of paper
# For white respondents:
mod.supp.W <- lm(V1_supp ~ V1_unfav + V1_fav, data=dataV1W_fair)
summary(mod.supp.W)
# Difference in support between favorable and equal = 0.09 (p<0.07) 
mod.supp.W2 <- lm(V1_supp ~ V1_equal + V1_fav, data=dataV1W_fair)
summary(mod.supp.W2)
# Difference in support between favorable and unfavorable = 0.65 (p<0.001) 

# For Black respondents:
mod.supp.B <- lm(V1_supp ~ V1_unfav + V1_fav, data=dataV1B_fair)
summary(mod.supp.B) 
# Difference in support between favorable and unfavorable = -0.01 (p<0.95) 

# Mediation analysis for white and Black respondents
# Figure 3 of main paper:  mediation analysis comparing favorable versus unfavorable agreements
# subset to favorable and unfavorable treatment

set.seed(43216)
dataV1fav.unfav <- subset(dataV1[dataV1$V1_unfav==1 | dataV1$V1_fav==1, ])
dataV1fav.unfav.W <- subset(dataV1fav.unfav[dataV1fav.unfav$white==1, ])
dataV1fav.unfav.B <- subset(dataV1fav.unfav[dataV1fav.unfav$white==0, ])

## Begin with mediation analysis for white respondents
## fairness of agreement as a mediator
fair3.med.W <- lm(V1_fair ~ V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav.W)
fair3.dv.W <- lm(V1_supp ~ V1_fair + V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav.W)
fair3.med2.W <- mediate(fair3.med.W, fair3.dv.W, treat="V1_fav", mediator="V1_fair", sims=1500)

# Extract quantities of interest for plot (note: plots can be made using the plot function, but they are displayed horizontally by default)
fair3.effects.W <- c(fair3.med2.W[1], fair3.med2.W[9], fair3.med2.W$tau.coef)
fair3.ci.low.W <- c(fair3.med2.W$d1.ci[1], fair3.med2.W$z0.ci[1], fair3.med2.W$tau.ci[1])
fair3.ci.high.W <- c(fair3.med2.W$d1.ci[2], fair3.med2.W$z0.ci[2], fair3.med2.W$tau.ci[2])

fair3.med.names <- c("ACME \n", "ADE \n ", "Total \n Effect")

## Now mediation analysis for Black respondents
## fairness of agreement as a mediator
set.seed(43216)
fair3.med.B <- lm(V1_fair ~ V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav.B)
fair3.dv.B <- lm(V1_supp ~ V1_fair + V1_fav + age + edu2 + income + party7 + women + NatAtt + ci, data=dataV1fav.unfav.B)
fair3.med2.B <- mediate(fair3.med.B, fair3.dv.B, treat="V1_fav", mediator="V1_fair", sims=1500)

# Extract quantities of interest for plot (note: plots can be made using the plot function, but they are displayed horizontally by default)
fair3.effects.B <- c(fair3.med2.B[1], fair3.med2.B[9], fair3.med2.B$tau.coef)
fair3.ci.low.B <- c(fair3.med2.B$d1.ci[1], fair3.med2.B$z0.ci[1], fair3.med2.B$tau.ci[1])
fair3.ci.high.B <- c(fair3.med2.B$d1.ci[2], fair3.med2.B$z0.ci[2], fair3.med2.B$tau.ci[2])

fair3.med.names <- c("ACME \n", "ADE \n ", "Total \n Effect")

# Combined plot of white and Black respondent mediation analysis
# Figure 3 of manuscript: export to PDF as 6x4
par(mar=c(3,5,1,2))
plot(1:3, fair3.effects.W , main = "", 
     ylab = "Average Support", xlab = "", cex.lab=2,
     ylim = c(-.35, .95), xlim = c(0.75, 3.25), pch = 22, type = "p", xaxt="n", cex=1)
mtext("", side=2, at=.2, cex=1.7, line=2.5)
abline(h=0, col="blue", lty=2)
for(i in 1:length(fair3.med.names)){
  lines(c(i,i), c(fair3.ci.high.W[i], fair3.ci.low.W[i]))
  mtext(fair3.med.names[i], side=1, at=i, 1.75, cex=1.15)
}
points(1.15:3.15, fair3.effects.B, pch=19, cex=1.5)
for(i in 1:length(fair3.effects.B)){
  lines(c(i+.15,i+.15), c(fair3.ci.high.B[i], fair3.ci.low.B[i]))
}
text(2.91, .93, labels = "Whites", cex=1.1)
text(2.9, .85, labels = "Blacks", cex=1.1)
points(2.5, .93, pch=22, cex=1.5)
points(2.5, .85, pch=19, cex=1.5)

# For quantities of interest cited in the text
summary(fair3.med2.W) # ACME 0.421 p<0.001, ADE 0.219, p<0.001, prop. mediated 0.656, p<0.001
summary(fair3.med2.B) # ACME 0.061 p<0.635, ADE 0.238, p<0.091, prop. mediated -0.031, p<0.968

# For effect cited in text on perception of fairness as equality as a mechanism
# Amongst those in the fairness experiment, test the effect of thinking fairness is about equality vs hard work

# The following line creates an indicator for those who thing of fairness as rewarding hard work
# Where this indicator equals zero, it means respondents think of fairness in terms of quality
dataV1_fair$reward.work <- ifelse(dataV1_fair$fair4==2, 1, 0)

mod.hard.work <- lm(V1_fair ~ V1_fav*reward.work, data=dataV1_fair[dataV1_fair$V1_equal!=1, ])
summary(mod.hard.work) # Interaction effect: 0.264 p<0.01


# National attachment as a moderator
# Replicate Table 1 Fairness Regressions controlling for nationalism

# Model for Table 2
mod.1.c <- lm(V1_fair ~ V1_equal*white +  V1_fav*white + NatAtt, data=dataV1_fair)
stargazer(mod.1.c)

# Nationalism interactions effects with treatment are reported in the appendix (Table 5)
# and are calculated in the appendix R code
 
# Test whether white liberals/Democrats adopt a principled fairness view, Table 3 of manuscript

# Ideology - lower values = liberal (scale is 1=extemely lib.... 5=extremely conservative)
# Party: 1=dem, 2=rep, 3=indep, 4=no pref

# Subset to white democrats
dataV1_W_Dem <- subset(dataV1W_fair[dataV1W_fair$party==1, ])
# subset to white republicans
dataV1_W_Rep <- subset(dataV1W_fair[dataV1W_fair$party==2, ])

mod.1.W.Dem <- lm(V1_fair ~ V1_fav, data=dataV1_W_Dem[dataV1_W_Dem$V1_equal!=1, ])

mod.1.W.Rep <- lm(V1_fair ~ V1_fav, data=dataV1_W_Rep[dataV1_W_Rep$V1_equal!=1, ])

# Subset to white liberals
dataV1_W_Lib <- subset(dataV1W_fair[dataV1W_fair$ideology<3, ])
# Subset to white conservatives
dataV1_W_Con <- subset(dataV1W_fair[dataV1W_fair$ideology>3, ])

mod.1.W.Lib <- lm(V1_fair ~ V1_fav, data=dataV1_W_Lib[dataV1_W_Lib$V1_equal!=1, ])

mod.1.W.Con <- lm(V1_fair ~ V1_fav, data=dataV1_W_Con[dataV1_W_Con$V1_equal!=1, ])

## Table 3 of manuscript
stargazer(mod.1.W.Dem, mod.1.W.Rep, mod.1.W.Lib, mod.1.W.Con)

